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<^ ; Abstract 



■ An effective Hamiltonian for the ferroelectric transition in PbTiOs is 
^ ' constructed from first-principles density-functional-theory total-energy and 

o ■ 

\G • linear-response calculations through the use of a localized, symmetrized ba- 

sis set of "lattice Wannier functions." Explicit parametrization of the polar 
lattice Wannier functions is used for subspace projection, addressing the is- 
sues of LO-TO splitting and coupling to the complementary subspace. In 
:'"! ' contrast with ferroelectric BaTiO^ and KNbOs, we find significant involve- 

■ ment of the Pb atom in the lattice instability. Monte Carlo simulations for 
this Hamiltonian show a first-order cubic-tetragonal transition at 660 K. Re- 
sulting temperature dependence of spontaneous polarization, c/a ratio and 
unit-cell volume near the transition are in good agreement with experiment. 
Comparison of Monte Carlo results with mean field theory analysis shows that 
both strain and fluctuations are necessary to produce the first order character 
of this transition. 
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I. INTRODUCTION 



Perovskite-structure oxides exhibit various types of lattice instabilities resulting from 
inherent structural frustration in the prototype cubic structureS, shown in Fig. 1(a). This 
class of materials includes a large number of ferroelectrics, with uniform polar distortions and 
accompanying lattice relaxation (eg. PbTiO-^, BaTiO-^, KNbO^), while cation substitution 
can result in dramatic changes in ground state distortion (eg. antiferroelectric in PbZrOs, 
antiferrodistortive in SrTiOs) and corresponding complexities in the mixed systems (eg. 
PhZri_xTixO^, Bai_xSrxTiO^). However, in nearly all examples, the amplitudes and en- 
ergies of the distortions are rather small, and cubic symmetry is restored at temperatures 
above a critical temperature Tc, typically a few hundred degrees Kelvin. 

For a better understanding of structural phase transitions in perovskite oxides, including 
chemical trends in the transition temperatures, the first-order vs second-order character of 
transitions, the relationship between local distortions and average crystallographic struc- 
ture, and the stability of intermediate temperature phases, first-principles calculations offer 
valuable access to microscopic information. With advances in algorithms and computational 
capabilities, the challenge of achieving the high accuracy necessary for studying these distor- 
tions has been largely met, and ground state distortions well reproduced for a wide range of 
perovskite-structure oxides^i. However, for ab initio molecular dynamics or Monte Carlo, 
the system sizes required for the study of finite-temperature structural transitions are still 
completely impractical. 

An alternative approach is to focus on a restricted subset of the degrees of freedom 
that is relevant to the transition and construct a simple effective Hamiltonian in this sub- 
space. The parameters in these models are chosen to reproduce the low energy surface 
of an individual material, and thus to reproduce its finite temperature behaviour in the 
vicinity of its transition. Comparison of these models gives a systematic understanding of 
similarities and differences in the microscopic structural energetics of different materials. 
From the dependence of calculated properties on effective Hamiltonian parameters, one can 
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also obtain a better understanding of the role of various microscopic couplings in producing 
the observed behaviour. Microscopic effective Hamiltonians for structural transitions in per- 
ovskites were first constructed using the concept of local modes, with empirically determined 
parameters.iS First-principles total energy calculations were used in the determination of a 
local mode effective Hamiltonian for the structural transition in GeTei and more recently 
for the structural transitions in PbTiO-^, BaTiO^ and SrTiO^. A systematic approach 
which generalizes and refines the local mode concept, allowing the efficient construction 
of an optimal effective Hamiltonian from first-principles total-energy and linear response 
techniques, has been developed based on the concept of lattice Wannier functions. This 
approach exploits symmetry properties of the system and is generally applicable to complex 
structural transitions involving several unstable modes including ones at the zone bound- 
ary. Information from additional first principles calculations allows for a systematic check 
on the validity of the truncation of the effective Hamiltonian, and, when needed, the ex- 
pansion of the subspace and refinements of its form. The resulting effective Hamiltonian is 
quantitatively realistic while retaining a simple and physically transparent form. 

In this paper, we present a detailed description of the first-principles investigation of 
PbTiO^, which exhibits a single first-order transition at 763 K from the cubic high- 
temperature phase to the ferroelectric tetragonal ground state, shown in Fig. 1(b). We 
construct an effective Hamiltonian for this structural phase transition from first principles 
using the lattice Wannier function method. In contrast with BaTiO^ and KNbO^ for 
which the uniform polar distortions in the low temperature phase consist of predominantly 
B-atom displacements, those in PbTiO^ are dominated by A-atom (Pb) displacements, which 
will be important in determining the effective Hamiltonian subspace. The effective Hamil- 
tonian also contains the coupling of these local polar distortions to strain. In Ref. ^, the 
tetragonal phase in PbTiO^ was found to be stabilized relative to the rhombohedral phase by 
the unit cell relaxation. In addition, we will find that strain plays a crucial role in producing 
the correct finite temperature transition behaviour. 

In Section IIA, we briefly review the method of lattice Wannier functions. In Section IIB 
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and lie, the first principles methods and results obtained for the lattice constant, elastic 
constants, phonon frequencies and the effective charges of PbTiO^ are presented. In Sections 
IIIA and IIIB, we describe the construction of the effective Hamiltonian, with particular 
attention to the treatment of LO-TO splitting and crossing of branches through explicit 
parameterization of the lattice Wannier functions. In Section IIIC, we describe properties 
of the ground state of the effective Hamiltonian determined from first principles. In Section 
IV, we present results of finite temperature analysis of Heff using mean field theory and 
Monte Carlo simulations. These results are discussed in Section V. 

II. METHOD 

A. Lattice Wannier Function method for the construction of i^e// 

In the lattice Wannier function method, the effective Hamiltonian is obtained as the re- 
sult of projection of the full lattice Hamiltonian (in the Born-Oppenheimer approximation) 
into a subspace of the full ionic displacement space. The effective Hamiltonian subspace 
is spanned by an orthonormal basis of "lattice Wannier functions:" symmetrized localized 
atomic displacement patterns taken with respect to a high-symmetry reference configura- 
tion. This basis defines a set of coordinates such that a given set of values of the coordinates 
corresponds directly to a particular pattern of atomic displacements. As a result of the 
symmetrized and localized nature of the basis, the Taylor expansion of the effective Hamil- 
tonian around the high-symmetry reference configuration (with all coordinate values equal 
to zero) has a simple form with relatively few parameters, which can be determined from 
first principles calculations using the correspondence to patterns of atomic displacements. 

We briefly review the procedure; further details can be found in Ref. [Tl|. Construction 
of the subspace begins with a Taylor expansion of the full lattice Hamiltonian to quadratic 
order. A subset of the eigenvectors of the quadratic Hamiltonian is selected for inclusion in 
the subspace. This subset must include the unstable modes which freeze in to produce the 
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low-temperature structure. In addition, to achieve a good description of the branches which 
emanate from the unstable modes, "end-points" of these branches at high symmetry k- 
points are included. The symmetry properties of the subspace are established by identifying 
symmetries of localized functions (Wyckoff position and site symmetry group irrep) which 
can build up the selected subset of modes. 



We follow the prescription in Ref. |TI] to obtain an explicit, though approximate, form 
of a lattice Wannier basis vector. This involves finding the symmetric coordination shells 
surrounding a representative Wyckoff site and identifying the independent displacement 
patterns of each shell that transform according to the given irreducible representation of the 
site symmetry group. The amplitudes of these displacement patterns completely specify an 
LWF. Because of the localized nature of LWFs, this infinite number of parameters can, to a 
good approximation, be reduced to a small finite number by neglecting the displacements of 
shells beyond a chosen range. At high symmetry points in the BZ, the modes built up with 
these parametrized LWFs are then fit to the corresponding normalized mode eigenvectors 
obtained from first principles. 

These basis functions completely specify the effective Hamiltonian subspace. In the 
ideal case, this subspace is completely decoupled at quadratic order from its complementary 
subspace and to a good approximation at higher order as well. This happens when the 
subspace consists of entire branches isolated in energy from the others and contains all 
the unstable modes. In most real systems, branches emanating from the unstable modes 
cross with branches in the complementary subspace. This leads to some degree of quadratic 
coupling which is unimportant if the crossing occurs far away from the unstable modes. 
If not, the subspace should be expanded to include these branches. In addition, in polar 
crystals, the electric field at g ^ can mix the LO modes differently from the corresponding 
TO modes. In such a case, the Wannier basis vector which reproduces a given TO branch 
will not reproduce any LO mode exactly. However, since LO modes are typically high in 
energy, this approximate description of the LO mode should be adequate for the description 
of the structural transition. 



The quadratic part of the lattice Hamiltonian has both kinetic and potential energy 
contributions. However, in the classical statistical mechanical treatment, the kinetic energy 
terms appear in Gaussian integrals in the partition function, factoring out to give a trivial 
contribution to the free energy. In this case, the eigenmodes used in the construction above 
could be of either the force constant matrix or the dynamical matrix. In PbTiO^, we have 
found that the difference in the resulting effective Hamiltonian subspace is rather small and 
both choices should give comparable results. In the construction described in section HI, 
we have used the eigenmodes of the force constant matrix. Eigenmodes of the dynamical 
matrix are strongly preferable only if the effective Hamiltonian is also to be used in classical 
dynamics or quantum mechanical simulations, since in that case the form of the kinetic 
energy is greatly simplified. 

In the final step, the lattice Hamiltonian is projected into this subspace to obtain the 
effective Hamiltonian. An explicit form H^ff is obtained by identifying a small number of 
physically important terms in a Taylor expansion in the lattice Wannier coordinates. The 
coefficients of these terms are parameters to be determined from first principles by fitting 
Heff to the results of selected total energy and linear response calculations, using the exphcit 
correspondence between the Wannier coordinates and the actual ionic displacements. To 
check the validity of the truncated form of the effective Hamiltonian, additional independent 
first principles calculations can be performed and compared with H^ff. 



B. First principles total energy calculations 

The first-principles calculations for PbTiO^ were performed using the ab initio pseu- 
dopotential method in the local density approximation (LDA) with the Perdew-Zunger 
parametrization of the Ceperley-Alder density functional.0 For Pb, the scalar-relativistic 
pseudopotentials from Ref. |1^ were used. The use of a plane wave basis set dictates the 
u.e of optimized p.eudopotentialJi for O and Ti to aAieve reasonable energy convergence 
and transferability. For O, the reference configuration 2s^2]9^ was used with pseudopotential 
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core radii Tc^s = r^^p = 1.5 a.u. Optimization was performed with gc,s = 7.0(-R?/)^/^ and 
<lc,p = 6-5{RyY^'^ and 4 and 3 Bessel functions for s and p orbitals respectively. For Ti, it 
is essential to treat the semi-core 3s and 3p electrons as valence electronsi'lll. The reference 
configuration 3s^3p^3cP was used with pseudopotential core radii r^s = r^p = 1.45 a.u. 
and Tc^rf = 1.5 a.u.. Optimization was performed with = 7.2(_R?/)^/^, gc,p = 7.0{RyY^^ 
and gc,d = 7.74(i?y)^/^, and 4 Bessel functions. An energy cutoff of 850 eV (corresponding 
to approximately 3600 plane waves for a 5-atom unit cell) was used to ensure convergence 
within 10mi??//atom. The self-consistent total energy calculations were performed using the 
program CASTEP 2.10 based on the stable and efficient preconditioned conjugate- gradients 
method^!. For the Brillouin zone integrations, k-point sampling was performed using the 
Monkhorst-Pack construction^ with 64 k-points in the full Brillouin zone. 

As reported in Ref. ^ and summarized here in Table (I), the lattice constant and elastic 
constants of PbTiO^ in the cubic perovskite structure obtained from the total energy calcula- 
tions for a range of unit cell volumes are in good agreement with previous calculationsi'i. In 
addition, in Fig. 4 of Ref. ^ we showed the calculated energies as a function of experimental 
soft mode amplitude, which compare favorably with previous LAPW calculation^. 



C. First principles DFT linear response 

The technique of DFT linear response is used to calculate the second derivatives of the 
total energy with respect to perturbation parameters through the self consistent calculation 
of the first order correction to the occupied Kohn-Sham wave function For example, 

Born effective charges, dielectric constant and dynamical matrices are the second derivatives 
of total energies and thus can be obtained with this technique. In this framework, the 
dielectric constant can be calculated avoiding cumbersome sums over unoccupied bands. 
Another significant advantage is that q ^ force constants can be computed with an effort 
similar to that of a single unit cell total energy calculation. 

Our implementation is a modification of CASTEP 2.1 based on the variational formula- 
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tion of DFT linear responseE3. All the linear response calculations reported here were done 
at the experimental lattice constanlH of 3.96883A with 64 Monkhorst-Pack k-points in the 
full Brillouin zone. The Vosko-Wilk-Nusair parametrization of the Ceperly-Alder density 
functional was used to permit the calculation of derivatives of exchange-correlation termsS. 
A 36 X 36 X 36 fourier transform grid is used for integration over a unit cell in real space. 
This real space grid breaks global translational invarianceil'ii, which manifests itself as small 
violations of the acoustic sum rule (the calculated frequencies of zone center acoustic modes 
are not exactly zero) and charge neutrality (the calculated change in polarization due to a 
rigid displacement of the crystal in any direction is not exactly zero). The acoustic sum 
rule was imposed by adding small corrections to the diagonal elements of the q = force 
constant matrix. Charge neutrality was imposed by adding the same small correction to the 
effective charges of all atoms. 

The Born effective charges are presented in Table (II) in very good agreement with 
previously calculated values^. The main features of interest are the anomalously large 
effective charges of Ti and O along the bond and the anisotropy of the oxygen charge. The 
calculated dielectric constant is 8.24, which can be compared with the experimental value of 
8.64 quoted in Ref. The data in Table (II) combined with the calculated force constants at 
q = give the frequencies of IR-active phonons presented in Table (III). Direct comparison 
of these results with the previous calculationJ^l is not possible because the calculations 
were performed at different lattice constants. This has an especially large impact on the 
unstable mode frequency, as confirmed by our calculations of coupling between this mode 
and homogeneous strain, to be described below. As can be seen in Table (IV), the unstable 
Fis mode has the largest mode effective charge, which should be associated with the largest 
LO-TO splitting. Since there are three polar zone center modes with the same symmetry, 
mixing leads to LO-mode eigenvectors different from TO-mode eigenvectors. Effects of this 
mixing can be quantified using the correlation matri> 



% =< e™|M|e^° >, (2.1) 



8 



given in table (IV), where Mmn = MmSmn is the mass matrix and are the IR-active mode 
eigenvectors. As expected, the unstable Fis TO mode has the strongest correlation with the 
highest LO mode. 

In Table V, we report selected phonon frequencies at other high symmetry k-points in 
the Brillouin zone, focusing in particular on the lowest energy phonons. While we also find 
unstable modes away from the F point, the unstable mode at F is clearly the dominant lattice 
instability in our calculations, consistent with the observed low-temperature structure. The 
eigenvectors of the lowest energy phonons and the corresponding force constant matrices 
will be used in the determination of the parameters in the effective Hamiltonian in the next 
section. 



III. CONSTRUCTION OF THE EFFECTIVE HAMILTONIAN FOR PBTI03 

A. Construction of the subspace 

The construction of the effective Hamiltonian subspace begins with consideration of 
the calculated force constant matrix eigenvalues and eigenvectors at F, R, M and X. The 
subspace has to include the unstable polar F15 mode which freezes in to produce the low tem- 
perature tetragonal structure. In addition, to achieve a good description of branches which 
emanate from this dominant unstable mode, the endpoints of these branches i?i5, Mg, M5, 
are included. As can be seen from Table (I) of Ref. |ll], the lattice Wannier functions which 
can build up this subset of modes transform like 3-dimensional vectors centered at Pb-sites. 
It should be noted that the lowest mode at R is actually R25, which corresponds to an 
oxygen octahedron rotation instability seen in many perovskite oxidesS. Since crossing of 
the lowest branch along (111) with that emanating from R25 occurs far from the relevant 
mode F15 and relatively higher in energy, we do not include it in the subspace. 

To include coupling of the relevant polar distortions (F15) to local distortions of the 
unit cell (inhomogeneous strain), we expand the subspace to include the acoustic modes by 
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choosing an additional set of lattice Wannier functions. Of the three possibilities (listed in 
Table (I) of Ref. |Tl|), Ti-centered 3- dimensional vectors are preferable to O^^x (1- dimensional 
vectors) and Ox,y (2-dimensional vectors), since this choice corresponds to the smallest 
subspace expansion and highest site symmetry group. Furthermore, unlike Ox^x-, the resulting 
6 ■ dimensional subspace does not include the highest energy modes. 

Next, we obtain an explicit form for the Pb-centered LWF. This involves finding the 
symmetric coordination shells surrounding a Pb site and identifying the independent dis- 
placement patterns of each shell that transform according to the vector representation of 
the site symmetry group Oh- For a given shell there can be more than one pattern of 
displacements with a given transformation property. To each such pattern corresponds an 
independent amplitude parameter. By neglecting the displacements of shells beyond 1st 
neighbour Ti and O shells and 2nd neighbour Pb shells, we obtain a total of 10 parameters. 
The first shell of Ti atoms has 2 independent displacement patterns, parametrized by ti and 
t2'i there are 1, 2, 2 parameters for the zeroth, first and second shells of Pb atoms respectively 
and 3 parameters for the first shell of oxygen atoms. These displacement patterns are shown 
in Fig. 2. 

To determine the numerical values of these parameters for PbTiO^, we build up the 
transverse modes Cq-^a at high symmetry k-points in the Brillouin zone, namely F, X, M, and 
R, from the parametrized LWF using 

eq,a = J2 exp{iq- Ri)Wi^a (3.1) 

where Ri is a direct lattice vector and Wi^a is an LWF centered at the Pb site in the ith unit 
cell, a being its Cartesian component. This specifies atomic displacements in these modes 
as linear functions of the parameters, to be fit to the normalized eigenvectors of the force 
constant matrix calculated from first principles. With the parameters listed above, we can 
reproduce the normalized eigenvectors of the modes Fis, i?i5 and exactly. The remaining 
free parameters, associated with Pb and Ti displacements, were used to fit to a normalized 
mode with maximum overlap with the lowest Mg (see Table (VI)). Numerical values of these 

10 



parameters, listed in Table (VII), clearly show that the magnitude of the parameters decays 
rapidly with shell-radius, confirming the assumption of LWF localization. Furthermore, by 
adding an additional shell of Pb (one parameter) and oxygen atoms (two parameters), we 
can reproduce all the transverse optical modes in the subspace. The parameter values for 
this refined LWF are given in Table (VII). The values of the parameters of the innermost 
shells do not change very much, and the values of the new parameters are very small. 

Another way of testing the approximate LWF is to see how well it reproduces other modes 
in the subspace. For example, in Table VIII, we show the comparison of the first principles 
X'q eigenvector with the mode constructed with the approximate LWF. The approximate 
mode has an overlap of 92% with the relevant mode, and if the approximate mode vector 
is normalized, the overlap becomes 99.96%, showing that the LWF describes the subspace 
very well. 

For the simplest treatment of inhomogeneous strain (the acoustic branches), an explicit 
expression of the Ti-centered LWF is not needed, since the goal is only to reproduce the 
long wavelength acoustic modes, whose dispersion is determined from the elastic constants. 
For a more refined treatment, an LWF could be parametrized as above and determined by 
fitting to the first principles eigenmodes ri5, i?25, M3, M5, Xi, X5. 

B. Determination of parameters in i?e// 

Using the symmetry properties of the lattice Wannier basis for the effective Hamiltonian 
subspace, we write an explicit expression for H^ff as a Taylor expansion in the lattice 
Wannier coordinates, invariant under the space group PmSm. {^i} and {ui} denote the 
Pb-centered and Ti-centered lattice Wannier coordinates respectively. Each of these three 
dimensional vector degrees of freedom transforms according to the Pis irrep of the point 
symmetry group Oh- Below, we organize the terms in the expansion of i^e// into those acting 
exclusively in the Pb-centered subspace and the Ti-centered subspace and those couphng 
the two. 
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In the Pb-centered subspace, we consider quadratic interactions up to third nearest 
neighbour with the most general form allowed by the space group symmetry: 

i 

+ E E Mii ■ Miid) ■ d) + ariii ■ iiid) - (6 • Miid) ■ d))] 

* d=nnl 

+ E E ■ Mi{d) ■ d)+bTl{ii ■ di){ii{d) ■ d,)+bT2{ii ■ d2){ii{d) ■ d2)] 

* d=nn2 

+ E E Mi^ ■ d)Ud) ■ d) + CT{i ■ Ud) - {i ■ d){Ud) ■ d))l (3.2) 

* d=nn3 

where ^i{d) denotes the LWF coordinate at a neighbour of site i in d direction. Beyond third 
neighbor we use a dipole-dipole form parametrized by the mode effective charge Z* and the 
electronic dielectric constant €00- 

• ^'^^^ ~ • ^^^^'^^^ (3.3) 
i J- Coo \d\^ 

An important simplifying approximation is that the onsite potential, depending on the value 

of at a single i, is the only set of terms including anharmonic interactions acting exclusively 

in the Pb-centered subspace. For simplicity, anharmonic terms are included only in the onsite 

potential with isotropic terms up to eighth order in and full cubic anisotropy at fourth 

order: 

+ C{^L + 4 + a + + E\m- (3.4) 

i 

In the Pb-centered subspace, the parameters to be determined from first principles are 
A,aL,aT,bL,bTi,bT2,CL,CT, B,C, D, E, Z . This determination relies on the explicit cor- 

— * — * 

respondence between the lattice Wannier coordinate {^j} and the ionic displacements {di} 
obtained in subsection III A. This correspondence allows us to relate the first principles total 
energies and the derivatives of total energies to various terms in i?e//- The parameters in the 
quadratic part of i^e// are linearly related to the force constant matrices obtained from den- 
sity functional linear response at high symmetry k-points in the BZ. In Table (IX) are given 
specific relations for modes at various k-points in the BZ including Fi^, X2,X'^, Mg, M2, R15 
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and the A3 modes at (lll)7r/(2a), (lll)7r/(4a). The parameter Z is determined from the 
eigenvector of unstable Fis and the Born effective charges. Obtaining eoo directly from DFT 
linear response and solving the system of linear equations yields values for all the parameters 
in the quadratic part of i?e//, listed in Table (IX). The resulting normal mode dispersion of 
ife// is shown in Fig. 3. For the LO modes at (lll)7r/(2a) and (lll)7r/(4a), the reasonable 
agreement between the calculated force constant matrix eigenvalue and i^e// is an indication 
of the validity of the truncation in the Taylor expansion. 

The parameters B.C. D.E appearing in the onsite anharmonic terms are determined 
from the total energies of uniformly distorted configurations (^^ = as shown in Fig. 
4. The minimum energy configuration has rhombohedral symmetry along (111)). The 
difference among the energies of uniform distortions with different symmetries ((100), (110), 
(111)) is a reflection of the cubic anisotropy, which is described quite well by the fourth 
order terms. The resulting parameters are listed in Table (X). 

To account for the effects of changes in lattice parameters at the structural phase tran- 
sition, we include the lowest order terms in the homogeneous strain and its coupling to the 
Pb-centered subspace: 

+5o(I]eaa)^|6r + 5l^(eaaI3Cj^a) +92 1^ e^/? (3-5) 
ce i a i a</3 i 

where is a component of the strain tensor, Cn, C12, C44 are the elastic constants, and the 
parameters go, 91,92 give the strength of coupling of strain with the local polar distortions 
^ia- All these parameters are determined from the single unit cell total-energy calculations 
for three independent types of unit cell distortions (isotropic, uniaxial and rhombohedral 
shear), with magnitudes of up to 2 to 4 % of the experimental lattice constants. The total 
energies of these strained unit cells with no local polar distortion, shown in Fig. 5, give 
the three elastic constants Cii,Ci2 and C44. For each of these unit-cell-strain types, we 
also compute the second derivative of energy with respect to uniform local polar distortions 

— * — * 

— ^, as shown in Fig. 5. These results yield the coupling parameters shown in Table (X). 
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Now we turn to the determination of the terms in Hejf acting in the Ti-centered subspace. 
Because this subspace contains the zone center acoustic modes, these terms must satisfy 
global translational and rotational invariance. This symmetry is built into the systematic 
expansion procedure given by KeatingS, in which invariant terms are built up from dot 
products of differences of the Wj's. If the expansion of the effective Hamiltonian in the 
Keating construction^ is truncated at quadratic order and three independent parameters, 
corresponding to the three elastic constants, the following terms are obtained: 

+ X! H [^L{ui ■ d){ui{d) ■ d) + driui ■ Ui{d) - {ui ■ d){ui{d) ■ d))] 

* d=nnl 

j^Y^ [biiui ■ d){ui{d) ■ d) + bTi{ui ■ di){ui{d^ (3.6) 

* d=nn2 

The relations of these parameters to the elastic constants are made by using the Keating 
expansion to evaluate the energy of homogeneously strained configurations. With these 
relations, A = C11+2C44, cll = -^Cu, dx = -^€44, hi = -hri = -|Ci2+^C'44 and 6t2 = 
0, these parameters can easily be obtained from first-principles calculations. Because there 
are no unstable modes in this subspace, there is no need to include higher-order interactions. 
In any case, within the local anharmonicity approximation, global translational invariance 
requires anharmonic terms to be zero at all orders. As mentioned in the previous section, 
there is no need for an explicit form of the Ti-centered LWF in this minimal treatment. 
For refinement of H^ff in this subspace, one could construct an explicit form and determine 
additional parameters in a manner analogous to that for the Pb-centered subspace. 

Finally, the simplest coupling between local polar distortions (Pb-centered subspace) and 
inhomogeneous strain (Ti-centered subspace) that satisfies the constraint of global transla- 
tional invariance and does not vanish in the limit A: — is the nearest-neighbor coupling 
linear in u and quadratic in ^, with both ^ variables taken on the same site: 

* d=±y±z 
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+ ^ E («^(^^ + f ^ + f ^1 - MR^ - f ^ + 1^1)) + C.J>.} 
d=±xzky 

+ ^ Y: MR^ + f ^ + - ^yiR^ - + 1^1)) + c-P-} (3-7) 

d=±y±z 

couples only to differences of the M/3's, which can be recognized as finite difference ap- 
proximations to the gradient, and thus as the local strain tensor (see also Ref. ^). The 
terms in H^jf coupling inhomogeneous strain and polar distortions are related in the long 
wavelength limit to the coupling between homogeneous strain and the polar Fis mode. Thus 
the three independent coupling parameters can be obtained from the corresponding homo- 
geneous strain coupling parameters: ho = {go + fi'i)/4, hi = go/4: and /i2 = 5'2/8. 



C. Examination of model energetics 

Having fully determined H^ff, we now explore the low energy surface of the model to 
confirm that it gives a correct ground state when compared with the real crystal. Since 
the anharmonic terms occur only in the Pb-subspace and are local (the anharmonicity 
is wavevector independent), it is easy to determine the ground state from the quadratic 
order terms. The lowest energy mode is obtained by freezing in the most unstable mode: 
Fis. We consider changes in energy as this mode is frozen in with polarization along the 
(001), (110) and (111) directions. In Fig. 4, it can be seen that the rhombohedral state 
((lll)-distortion) has the lowest energy. If the unit cell is allowed to relax as the mode is 
frozen in, by minimizing over the homogeneous strain, we find an overall increase in distortion 
energy, with the lowest energy state being of tetragonal symmetry ((OOl)-distortion) as can 
be seen in Fig. 6. This is consistent both with previous first principles calculationsi and 
experimental resultsii. 

For the lowest energy tetragonal configuration, we obtain a value for the spontaneous 
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polarization from the mode effective charge of Ps = 0.87C/m^, which is in the range of values 
{Ps = 50 to lOOC/m^) reported from experimental. From the values of homogeneous strain 
in this ground state, we obtain a c/a ratio of 1.08, to be compared with the experimental 
value of l.Oeii. Using the explicit form for the LWF, the atomic displacements in the model 
ground state can be obtained. We find that the oxygen octahedra are almost undistorted and 
the relative displacement of the Pb atoms is 0.35 A(to be compared with the experimental 
value of 0.47 Ai). It should be noted that this is not the true LDA ground state, since 
previous total energy calculations! for the experimental distortions showed the latter is 
slightly lower in energy. This means that there is a higher order coupling of the unstable Fis 
mode to an additional polar F15 mode not included in the effective Hamiltonian subspace. 
Considering that the atomic displacements in PbTiO^ are relatively large, the presence of 
such anharmonicity is not surprising. However, as discussed above, the ground state of the 
model is very similar to the experimental ground state and the small loss of accuracy is more 
than compensated for by the gain in simplicity. 



IV. FINITE-TEMPERATURE BEHAVIOUR 

The effective Hamiltonian is constructed to show the same finite temperature critical 
behaviour as the full lattice Hamiltonian in a statistical mechanical analysis. While the form 
of Heff is somewhat too complicated for the application of analytical methods beyond mean 
field theory, it is quite suitable for Monte Carlo simulations, since the changes in energy for 
changing system configurations can be readily evaluated. Monte Carlo simulations are used 
in the detailed analysis of H^jf to obtain T-dependence of a variety of structural properties 
near the transition, while our mean field theory analysis is limited to the estimation of Tc 
and the identification of the order of the transition and symmetry of the phases. Comparison 
of the mean field results with those of Monte Carlo simulations allows us to study the effects 
of fluctuations. 
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A. Mean field theory 



Variational mean field theory for the class of models with variable length vector degrees 
of freedom and strain coupling was developed in Ref. ^ In this approach, the homogeneous 
strain and uniform polarization (P) are identified as the order parameters for the transition. 
P is directly related to the average value of uniform local distortion through the mode 
effective charge Z* and the unit cell volume ^ceih 

P = r. <e>Me//. (4.1) 

In the high temperature (cubic perovskite) phase, the uniform polarization is zero and 
the strain tensor has full cubic symmetry: Cxx = ^yy = Gzz- We used the variational 
formulation of mean field theory, which involves constructing a trial density matrix as a 
product of single site density matrices and minimizing the resulting free energy functional 
with respect to the variational parameters in the trial density matrix. We minimized the 
trial free energy with respect to variational parameters corresponding to cubic, tetragonal 
and rhombohedral symmetries to determine the stable phase at various temperatures. The 
system is stable in the cubic phase above the transition temperature T^. = llOOi^' and in 
the tetragonal phase below T^, but within the accuracy of our calculation, the transition is 
second order. Switching off the coupling to homogeneous strain resulted in a second order 
cubic-rhombohedral transition at a significantly lower temperature of 900 K. 

B. Monte Carlo simulations 

Classical Monte Carlo simulationsH were performed using the Metropolis algorithm for 
finite size systems of L x Lx L unit cells and periodic boundary conditions. A configuration 
of the system is specified by two sets of three dimensional vectors {C,i} and {uj} placed on 
interpenetrating simple cubic lattices of size L x L x L. We generated a trial configuration 
by updating a single vector to a randomly chosen vector inside a cubic box centered at the 
current value of the vector. The size of this box is chosen to yield a reasonable acceptance 
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ratio (> 0.1), and is roughly 0.25ao near Tc. With the change in a single vector, the change in 
energy associated with short range terms (quadratic interactions up to third neighbour, the 
onsite potential, coupling to strain and third order coupling between the two subspaces) is 
easy to calculate. Because of their long range, computation of dipolar intersite interactions 
is relatively costly, limiting the size of our simulations to L < 12. The 3x3 matrix of 
dipolar intersite coupling for each pair of spins was calculated using the Ewald summation 
technique for each value of L. Changes in the quadratic intersite interaction due to changes 
in strain are neglected in this model. One Monte Carlo sweep (mcs) involves one update of 
the ^i's (in typewriter mode) followed by one update of the Wj's (in typewriter mode), and 
20 updates of the 6 components of the strain tensor. 

Preliminary Monte Carlo simulations performed with 25,000 to 50,000 sweeps showed 
dependence on the initial configuration of simulations at temperatures in the vicinity of the 
transition for L > 5. For L = 5, which is small enough for ergodic sampling in a run of this 
length, the energy histogram shows two clearly separated peaks. This behaviour is typical 
of a first order transition^ At larger system sizes, due to the exponentially increasing 
free energy barrier between the regions of configuration space corresponding to low and 
high temperature phases, only one of the two peaks in the energy histogram is sampled, 
depending on the choice of initial configuration. An accurate determination of Tc requires 
knowledge of the relative free energies of the high and low temperature phases as a function 
of temperature. Recently developed methods to extract for first order transitions include 
multicanonical algorithms!^ and jump-walking algorithmslll. In our applications of these 
methods to PbTiO^, we found that these approaches require very long (10^ mcs) simulations, 
and therefore are rather impractical. However in the present case, the uncertainty in Tc 
obtained from the range of hysteresis is small compared to the LDA and other errors in our 
analysis and therefore high accuracy determination of is not necessary. The calculation 
of the physical properties of the high and low temperature phases at temperatures inside 
the range of hysteresis is carried out with an appropriate choice of the initial configuration. 

In Fig. 7, we show the bounds on Tc for L = 5 through 11, obtained by monitoring the 
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sensitivity of the average structural parameters to the choice of initial state: T> is the lowest 
temperature at which the system averages are characteristic of the cubic state, starting with 
an inital ground state tetragonal configuration, while T< is the highest temperature at which 
a starting cubic configuration results in a tetragonal state. A value of Tc = 660K, obtained 
from averaging the bounds at the largest system size, is in very good agreement with the 
experimental transition temperature 763K. 

To detect the symmetry of the low temperature phase, we calculated the averages of 
largest, smallest and intermediate absolute values of the cartesian components of ^ = 
These averages for the 7x7x7 system as a function of temperature are shown 
in Fig. 8. Near Tc, the largest component jumps to a finite value, while the other two 
components remain close to zero, indicating tetragonal symmetry of the low temperature 
phase. As shown in Fig. 9, this uniform tetragonal polar distortion is accompanied by a 
tetragonal strain c/a ^ 1. This quantity also shows a marked jump near Tc. Finally, from 
the average homogeneous strain we obtained the average volume of the system as a function 
of temperature, as shown in Fig. 10. The negative thermal expansion in the simulations 
just below Tc is also seen experimentallyii. 

The latent heat of a first order transition is given by the difference in energies at which 
the two peaks appear in the energy histogram in the simulations at Tc. To determine this 
energy difference, we performed two simulations for L = 9 at the midpoint of the hysteresis 
range Tc = 675K, one starting with a tetragonal configuration and the other starting with a 
cubic configuration. The difference in the positions of the peaks in the energy histograms for 
these two simulations yields an estimate of latent heat of 3400 J/mol, in rough agreement 
with the measured value of 4800 J/moleEl. It should be noted that these values are much 
larger than 209 J/mol latent heat of the cubic tetragonal transition of BaTiO^. 

Information about the local distortions in the high temperature nonpolar phase just 
above Tc, can be obtained from the single spin distribution < >. For all L, we find the 
distribution to be very close to Gaussian. The rather broad width 0.04ao) shows that 
there are significant local distortions. 
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For the system sizes used in the simulations, the couphng to the inhomogeneous strain 
appears to be relatively unimportant. If the coupling is set to zero the changes in the calcu- 
lated Tc and other transition properties are negligible. However, for larger scale simulations 
involving multiple domains the effects should be significant. 



V. DISCUSSION 

The first principles effective Hamiltonian constructed in the previous section provides 
a quantitative microscopic description of the structural energetics of PbTiO^ relevant to 
the paraelectric-ferroelectric phase transition. This model can be used to investigate the 
connection of specific features of the Hamiltonian to the observed behaviour in the vicinity 
of the transition. In addition, it is possible to connect these features to aspects of the 
chemistry of PbTiO^ and related compounds. 

One important feature of the Hamiltonian is that the TO branches are unstable through- 
out most of the BZ (Fig. 3). So, although Fis is the dominant instability, finite wavelength 
fluctuations have relatively low energy. This may account for the breadth of the single site 
distribution, and can be expected also to be reflected in the short range order. The unstable 
branch along the (111) direction in PbTiO^ is quite fiat when compared with SrTiO^ 
and KNboM- 

In comparison with KNbO-i and ferroelectric BaTiO^, in which the polar 
unstable modes have a strong B-componentlili, the instabilities in PbTiO^ are dominated 
by large Pb-displacements. From symmetry argument^ll, the Pb displacements couple with 
oxygen displacements at F, X, R and M leading to the low energy of the modes at those 
points. While the same argument applies to Ti displacements at T,X and M, they do not 
couple with any other mode at R point. Therefore, the energy of the mode at R is high and 
the dispersion along F to R is large. The special role Pb plays in the instabilities of PbTiO^ 
and PbZrO^ in contrast with the A-atoms in non-P6 perovskite compounds has its origin 
in the strong hybridization of Pb with oxygen atomsi'0. 

To understand the consequences of the couphng of the polar subspace to the strain at 
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finite temperature, we performed Monte Carlo simulations for H^ff with this coupling set to 
zero. This corresponds to a constant volume phase transition with the unit cell constrained 
to be cubic. In this case, we find a second order phase transition at 400 K directly to the 
rhombohedral phase. Thus, the coupling of local polar distortions to strain is responsible 
for both the stability of the tetragonal phase relative to the rhombohedral one and the first 
order character of the transition at finite temperature. As discussed in the previous section, 
mean field theory shows a second order transition both with and without strain coupling, 
implying that both fluctuations and strain coupling are required for producing the first order 
transition. Comparing the transition temperatures obtained in Monte Carlo and mean field 
theory with and without strain coupling, we find that while fluctuations suppress T^, the 
coupling to strain enhances the stability of the ferroelectric phase. 

In comparison with related ferroelectric compounds, the transition in PhTiO^ has a much 
stronger first order character, reflected in its large latent heat. While the strain coupling is 
responsible for the first order character^, anharmonicity in the lattice plays an important 
role in the magnitude of its discontinuityil. The minimum energy uniform polar distortions 
in PhTiO^ are much larger than those in related compounds indicating a large contribution 
from anharmonicity in the low-energy surface. The relation of these features to the chemistry 
of A or B atoms was discussed using ionic radii in Ref. 

There are two main sources of error in the work presented in this paper. One of these 
is the LDA used in the exchange correlation functional. Equilibrium lattice constants are 
typically underestimated in the LDA calculations. This can strongly affect the study of 
structural phase transition, since the lattice instabilities are very sensitive to the lattice 
parameters. In the present work these errors were partially eliminated by expanding the 
effective Hamiltonian around the experimental lattice constant near and dropping the 
term linear in strain. The other source of error is the truncation of the effective Hamilto- 
nian subspace. In the LWF approach, this subspace is decoupled at quadratic order from 
its complementary subspace to a good approximation. However there can be anharmonic 
coupling between the two subspaces. In the case of PbTiOs, there is a small higher order 
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coupling of Fis to modes not included in the subspace, which affects the energies of large 
distortions. Since these large distortions are mainly important at low temperatures, we 
expect this coupling to have a small effect on Tc- 

VI. CONCLUSION 

In conclusion, we have applied the method of lattice Wannier functions to construct an 
effective Hamiltonian for the ferroelectric phase transition in PbTiO^ completely from first 
principles. Monte Carlo simulations for this Hamiltonian yield a first order cubic-tetragonal 
transition at 660 K and a description of the system near the transition in good agreement 
with experiment. The strong involvement of Pb atom in the lattice instability as well as 
anharmonicity and the coupling of polar distortions with homogeneous deformations of the 
lattice are found to be very important in producing the transition behaviour characteristic 
of PbTiOa. 
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TABLES 

TABLE L Cubic perovskite lattice and elastic constants calculated from various first principles 
calculations. Elastic constants are given in eV/cell. 



This work Ref. | Ref. | 



ao (A) 3.883 3.889 3.889 

B (GPa) 203 215 209 

Cii 117. - 123. 

Ci2 51.6 - 53.6 

C44 137. - 148. 



TABLE IL Effective charges calculated from first principles linear response and compared with 
the results of the geometric phase approach (Ref. p6[). 



z. 



ol 



z. 



o2 



This work 



Ref. 26 



3.87 
3.90 



7.07 
7.06 



-5.71 
-5.83 



-2.51 
-2.56 



TABLE in. IR active optical phonon frequencies (cm^^) at F obtained using linear response 
at the experimental volume. They are compared with the results of the frozen phonon calculations 
performed at the LDA volume with ultrasoft pseudopotentials0. 









TOl 


T02 


T03 


LOl 


L02 


L03 


Pres( 
Ref. 


3nt 
26 


work 


182 I 
144 I 


63 
121 


447 
497 


47 
104 


418 
410 


610 
673 
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TABLE IV. LO-TO splitting: mode effective charges and correlation matrix. 



T LOl L02 LOS 



TOl 


9.45 




0.224 


0.466 


0.855 


T02 


2.56 




0.974 


0.116 


0.192 


T03 


1.53 




0.010 


0.876 


0.481 


TABLE V. 


Selected phonon 


frequencies (cm 


at high symmetry k-points calculated using 


DFT linear response. Symmetry labels follow the convention of Ref. 11. 




k-point 




phonon 






frequencies 


X 










30.6 I, 264 












93.1, 647 


M 










35.1 I, 400, 201 












16.4 












145 I 












15.5, 339 


R 




R25 






367 






R'u 






370 






R'2 






746 






Ai 






8.78, 249, 421, 696 


(111)^ 




A2 






148 






A3 




58.2 I, 82.9, 230, 301, 430 
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TABLE VI. Determination of LWF parameters. Linear combinations of these parameters for 



the modes at high symmetry k-points and the corresponding 


; components of the normahzed eigen- 


vectors of the force constant matrix. 




Mode 


rmnliiTiri.tioTi of thf naTRiriPtprs 

R_/.i.j.j.cA> u.i.vyj.J- vy-L u j-lv^ i_/ (a>x ctiJ-J-iv^ u v^x \j 


component of the eigenvector 


J- 15 


/r-t _l_ Am _l_ Om _L 1 On . 

Pi ~r ^P2 ~r ^Ps ~r l-^PA. 


0.5560 






0.5375 




■^'^ J. 


-0.3414 




402 


-0.4109 


Rib 


Pi - 4p2 - 2p3 + 12p4 


0.8981 




403 - 806 


-0.3110 




Pi - 4p2 + 2p3 - 4p4 


1.0000 




Pi - 2p3 - 4p4 


0.9010 




8*2 


0.3024 
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TABLE VII. Values of the LWF parameters determined from first principles. The parameters 
of the approximate LWF described in the text are given in the second column. Parameters for 
the refined LWF are obtained by fitting to all the TO modes at F, R, X and M, with additional 
parameters associated with third neighbour shell of Pb atoms and second neighbour shell of oxygen 

atoms. 



Parameter 


Approx. LWF 


Refined LWF 


Pi 


0.839 


0.829 


P2 


-0.037 


-0.049 


P3 


-0.012 


0.014 


P4 


-0.009 


-0.019 


P5 


0.0 


0.017 


Ol 


-0.085 


-0.086 


02 


-0.102 


-0.103 


03 


-0.077 


-0.087 


04 


0. 


.00045 


05 


0. 


-.0045 


h 


0.067 


0.067 


t2 


0.038 


0.037 



TABLE VIII. Comparison of eigenvectors. Mode vector (first row) built up using the 



approximate LWF is compared with the eij 


jenvector (second row) 


of the force constant matrix at 


X. 








Pb component 


component 


Mode in the subspace 


0.853 


-0.341 


Eigonvoctor from LR 


0.937 


-0.349 
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TABLE IX. Determination of coefficients in the quadratic part of Hgff. Linear combinations 
of these coefficients for the modes in the Hgff subspace at high symmetry k-points are equated to 
the corresponding eigenvalues of the projected force constant matrix. 





IV/Tr^Hp PicrpTivPilnp at nf fViP pfFppl'ivp TTflmilf r^niPin 


A/flliip ■frmn TiR 




7 — "Z*^ /(= 
Z — Zj /too 


19 18 


J- 15 




1 QflS 


X' 


A — 9/1 r 4- 4/7^7^ — d(hT 4- /iT-i "1 4- dhrrn — fi(rr 4- Pr-T-'l 4- 9 9^1 ^QQr /9 


fi 4fi7 


X' 
^5 








^ - 2aL - 46t2 + 8(cL + 2cr)/3 + 0.6165696z/2 


-0.360 




^ + 2aL - 4aT - 4(6l + 6ti) + 46t2 + 8(cl + 2ct)/3 - 1.23314z/2 


0.103 


Rib 


A - 2aL - 4aT + 4(6l + 6ti) + 46t2 - 8(cl + 2ct)/3 


0.076 


(111)^ 


A - {-2{bL - bxi) + 0.41635523z/2.0) 


-0.568 




A + V2(aL + 2aT) + 2(6l + bri + 6t2) + 0.942809(cl + 2ct) 






-((-^'L + bri) - 0.942809(cl - ct) + 0.7953677^/2) 


-1.750 



TABLE X. Parameters in the effective Hamiltonian (units of eV per unit cell) . 



A 


18.43 




39.27 


Cn 


117.9 


B 


2.629x10^ 




-10.67 


Cl2 


51.50 


C 


4.277x10^ 


hL 


4.882 


C44 


137.2 


D 


-1.658x10^ 




-1.391 


50 


-107.7 


E 


9.630x10^ 


bT2 


-0.1434 


91 


-790.3 


-7^*2 1 

Z /eoo 


12.18 


CL 


-3.389 


92 


-357.09 






CT 


0.7104 


f 


4.48 



30 



FIGURES 

FIG. 1. (a) Unit cell of the cubic perovskite compounds ABOs (b) Low temperature crystal 
structure of PbTiO^. Displacements of the atoms indicated by arrows form the polar distortions 
of the cubic unit cell. 

FIG. 2. z component of the vector-like Pb-centered lattice Wannier functions. P6, Ti and O 
atoms are represented by solid squares, empty squares and circles respectively. Parameters labeling 
the displacement patterns correspond to the length of the displacements (arrows) of atoms for the 
unit value of the LWF coordinate. 

FIG. 3. Normal mode dispersion of Heff. Solid circles are the first principles mode eigenvalues 
used in the fitting. Open circles are the first principles mode eigenvalues not used in fitting the 
Heff, which test the validity of the truncated form of the effective Hamiltonian. 

FIG. 4. Total energies for uniformly distorted configurations (^j = ^) along directions 
(001), (110) and (111). Solid lines are the fit obtained with the parameters B,C,D and E in 

Heff. 

FIG. 5. Energetics of the homogeneous strain ( (a) isotropic, (b) uniaxial and (c) shear) and 
its coupling to the uniform polar distortions. Circles are the total energies for the strained unit cell 
configurations with no polar distortions. Solid lines going through the circles are the fits obtained 
with the elastic constants Cii,Ci2 and C44. Squares correspond to the second derivative of the 
total energies with respect to uniform polar distortions for the strained unit cells. Solid lines going 
through the squares are the fits obtained with the coupling parameters go,9i and g2- 

FIG. 6. Model energetics of the uniform polar distortions along (100), (110) and (111). Dotted 
lines correspond to the polar distortions with the unstrained cubic unit cell, and solid lines to the 
distortions with unit cell allowed to relax with respect to homogeneous strain. 

FIG. 7. Bounds on the transition temperature Tc as a function of system size used in the 
simulations. 
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FIG. 8. Averages of largest, smallest and intermediate absolute values of the cartesian com- 
ponents of the order parameter ^ = ^ ^ function of temperature, obtained from Monte 
Carlo simulations (125000 mcs) of a 7 x 7 x 7 system. 

FIG. 9. Average tetragonal strain parameter ^ as a function of temperature, obtained from 
Monte Carlo simulations (125000 mcs) of a 7 x 7 x 7 system. 

FIG. 10. Average unit cell volume as a function of temperature, obtained from Monte Carlo 
simulations (125000 mcs) of a 7 x 7 x 7 system. 
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